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Abstract 



We review various algorithms for finding apparent horizons in 3 + 1 numerical 
relativity. We then focus on one particular algorithm, in which we pose the 
apparent horizon equation H = \I + KijV^n? — ii' = as a nonlinear 
, elliptic (boundary-value) PDE on angular-coordinate space for the horizon 

\ shape function r = h{9,(j)), finite difference this PDE, and use Newton's 



method or a variant to solve the finite difference equations. 



OO ' We describe a method for computing the Jacobian matrix of the finite 



differenced H{h) function H(h) by symbolically differentiating the finite dif- 
, ference equations, giving the Jacobian elements directly in terms of the finite 

Q I difference molecule coefficients used in computing H(h). Assuming the finite 

differencing scheme commutes with linearization, we show how the Jacobian 
elements may be computed by first linearizing the continuum H{h) equa- 
tions, then finite differencing the linearized continuum equations. (This is 
essentially just the "Jacobian part" of the Newton-Kantorovich method for 
solving nonlinear PDEs.) We tabulate the resulting Jacobian coefficients for 
' a number of different H(h) and Jacobian computation schemes. 

We find this symbolic differentiation method of computing the H(h) Jaco- 
bian to be much more efficient than the usual numerical-perturbation method, 
and also much easier to implement than is commonly thought. 

When solving the discrete H(h) = equations, we find that Newton's 
method generally shows robust convergence. However, we find that it has a 
small (poor) radius of convergence if the initial guess for the horizon posi- 
tion contains significant high-spatial-frequency error components, i.e. angular 
Fourier components varying as (say) cosmO with m ^ 8. (Such components 
occur naturally if spacetime contains significant amounts of high-frequency 
gravitational radiation.) We show that this poor convergence behavior is not 
an artifact of insufficient resolution in the finite difference grid; rather, it 
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appears to be caused by a strong nonlinearity in the continuum H{h) func- 
tion for high-spatial-frequency error components in h. We find that a simple 
"line search" modification of Newton's method roughly doubles the horizon 
finder's radius of convergence, but both the unmodified and modified meth- 
ods' radia of convergence still fall rapidly with increasing spatial frequency, 
approximately as \jrr?l'^ . Further research is needed to explore more robust 
numerical algorithms for solving the H(h) = equations. 

Provided it converges, the Newton's-method algorithm for horizon find- 
ing is potentially very accurate, in practice limited only by the accuracy of 
the H(h) finite differencing scheme. Using 4th order finite differencing, we 
demonstrate that the error in the numerically computed horizon position 
shows the expected O((A0)^) scaling with grid resolution A0, and is typi- 
cally -10-^(10"*^) for a grid resolution of M = ^(f^). 

Finally, we briefiy discuss the global problem of finding or recognizing 
the outermost apparent horizon in a slice. We argue that this is an impor- 
tant problem, and that no reliable algorithms currently exist for it except in 
spherical symmetry. 

04.25.Dm, 02.70.Bf, 02.60.Cb, 02.60.Lj 
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Typeset using REVT^ 



I. INTRODUCTION 



In 3 + 1 numerical relativity, one often wishes to locate the black hole(s) in a (spacelike) 
slice. As discussed by Refs. [p3| , p4| , a black hole is rigorously defined in terms of its event 
horizon, the boundary of future null infinity's causal past. Although the event horizon has, 
in the words of Hawking and Ellis (Ref. ||3^), "a number of nice properties", it's defined in 
an inherently acausal manner: it can only be determined if the entire future development of 
the slice is known. (As discussed by Refs. 0,0], in practice the event horizon may be located 
to good accuracy given only the usual numerically generated approximate development to a 
nearly stationary state, but the fundamental acausality remains.) 

In contrast, an apparent horizon, also known as a marginally outer trapped surface, 
is defined (Refs. ||5^J5^) locally in time, within a single slice, as a closed 2-surface whose 



outgoing null geodesies have zero expansion. An apparent horizon is slicing-dependent: if we 
define a "world tube" by taking the union of the apparent horizon(s) in each slice of a slicing, 
this world tube will vary from one slicing to another. In a stationary spacetime event and 
apparent horizons coincide, although this generally isn't the case in dynamic spacetimes. 
However, given certain technical assumptions, the existence of an apparent horizon in a slice 
implies the existence of an event horizon, and thus by definition a black hole, containing 
the apparent horizon. (Unfortunately, the converse doesn't always hold. Notably, Wald and 
Iyer (Ref. [^) have constructed a family of angularly anisotropic slices in Schwarzschild 



spacetime which approach arbitrarily close to r = yet contain no apparent horizons.) 

There is thus considerable interest in numerical algorithms to find apparent horizons in 
numerically computed slices, both as diagnostic tools for locating black holes and study- 
ing their behavior (see, for example, Refs. 0,0), and for use "on the fiy" during numeri- 
cal evolutions to help in choosing the coordinates and "steering" the numerical evolution 
(Refs. p7| , p2| , p8| ,p| ) . This latter context makes particularly strong demands on a horizon- 
finding algorithm: Because the computed horizon position is used in the coordinate condi- 
tions, the horizon must be located quite accurately to ensure that spurious finite difference 
instabilities don't develop in the time evolution. Furthermore, the horizon must be re- 
located at each time step of the evolution, so the horizon-finding algorithm should be as 
efficient as possible. Finally, when evolving multiple-black-hole spacetimes in this manner 
it's desirable to have a means of detecting the appearance of a new outermost apparent 
horizon surrounding two black holes which are about to merge. We discuss this last problem 
further in section Xl. 

In this paper we give a detailed discussion of the "Newton's method" apparent-horizon- 
finding algorithm. This algorithm poses the apparent horizon equation as a nonlinear elliptic 
(boundary-value) PDE on angular-coordinate space for the horizon shape function r = 
h{6,(f)), finite differences this PDE, and uses some variant of Newton's method to solve 
the resulting set of simultaneous nonlinear algebraic equations for the values of h at the 
angular-coordinate grid points. This algorithm is suitable for both axisymmetric and fully- 
general spacetimes, and we discuss both cases. As explained in section y, we assume a 
locally polar spherical topology for the coordinates and finite differencing, though we make 
no assumptions about the basis used in taking tensor components. 
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II. NOTATION 



Our notation generally follows that of Misner, Thorne, and Wheeler (Ref. f^), with 
G = c = 1 units and a (—,+,+,+) spacetime metric signature. We assume the usual 
Einstein summation convention for repeated indices regardless of their tensor character, and 



we use the Penrose abstract-index notation, as described by (for example) Ref. [^. We use 



the standard 3 + 1 formalism of Arnowitt, Deser, and Misner (Ref. 0) (see Refs. [^,^] for 
recent reviews). 

We assume that a specific spacetime and 3 + 1 (spacelike) slice are given, and all our 
discussions take place within this slice. We use the term "horizon" to refer to the (an) 
apparent horizon in this slice. We often refer to various sets in the slice as being 1, 2, or 
3-dimensional, meaning the number of spatial dimensions - the time coordinate is never 
included in the dimensionality count. For example, we refer to the horizon itself as 2- 
dimensional. 

We assume that the spatial coordinates = (r, 6, (p) are such that in some neighborhood 
of the horizon, surfaces of constant r are topologically nested 2-spheres with r increasing 
outward, and we refer to r as a "radial" coordinate and 6 and as "angular" coordinates. For 
pedagogical convenience (only), we take 6 and to be the usual polar spherical coordinates, 
so that if spacetime is axisymmetric (spherically symmetric), is {6 and are) the symmetry 
coordinate(s). However, we make no assumptions about the detailed form of the coordinates, 
i.e. we allow all components of the 3-metric to be nonzero. 

We emphasize that although our assumptions about the local topology of r are funda- 
mental, our assumptions about the angular coordinates are for pedagogical convenience only, 
and could easily be eliminated. In particular, all our discussions carry over unchanged to 
multiple black hole spacetimes, using (for example) either Cadez conformal- mapping equipo- 
tential coordinates (Ref. |T6|) or multiple-coordinate-patch coordinate systems (Ref. p6|). 

We use ijkl for spatial (3-) indices, and uvwxy for indices ranging over the angular coor- 
dinates only. Qij denotes the 3-metric in the slice, g its determinant, and Vj the associated 
3-covariant derivative operator. Kij denotes the 3-extrinsic curvature of the slice, and K its 
trace. 

We use A to denote the 2-dimensional space of angular coordinates {0, 0). We sometimes 
need to distinguish between field variables defined on A or on the (2-dimensional) horizon, 
and field variables defined on a 3-dimensional neighborhood M of the horizon. This distinc- 
tion is often clear from context, but where ambiguity might arise we use prefixes and 
respectively, as in ^"^^H and '^^'>H. 

We use italic Roman letters h, etc., to denote continuum coordinates, functions, 
differential operators, and other quantities. We use sans serif Roman letters H, h, etc., to 
denote grid functions, and small capital Roman indices i, j, and k to index grid points. We 
use subscript grid-point indices to denote the evaluation of a continuum or grid function at 



a particular grid point, as in Hi or Hj. We use J P(Q) 



to denote the Jacobian matrix of 



product of two 
to denote the 



the grid function P = P(Q), as defined by Eq. (|T^, and ■ to denote the 
such Jacobians or that of a Jacobian and a grid function. We use J P{Q) 
linearization of the differential operator P = P{Q) about the point Q. 

We use M as a generic finite difference molecule and m as a generic index for molecule 
coefficients. We write m G M to mean that M has a nonzero coefficient at position m. 
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Temporarily taking (m) to denote some particular coordinate component of m, we refer to 
maxjvieM 1(^)1 as the "radius" of M, and to the number of distinct (m) values with m G M as 
the "diameter" or "number of points" of M. (For example, the usual symmetric 2nd order 
3-point molecules for 1st and 2nd derivatives both have radius 1 and diameter 3.) We often 
refer to a molecule as itself being a discrete operator, the actual application to a grid function 
being implicit. 

Given a grid function / and a set of points {xk} in its domain, we use interp(/(x), x = a) 
to mean an interpolation of the values f{xk) to the point x = a and interp'(/(a;), x = a) 
to mean the derivative of the same interpolant at this point. More precisely, taking I to be 
a smooth interpolating function (typically a Lagrange polynomial) such that /(x^) = f{xk) 
for each fc, interp(/(x), x = a) denotes /(a) and interp'(/(x), x = a) denotes (9//(9x) 

x=a 

III. THE APPARENT HORIZON EQUATION 

As discussed by (for example) Ref. [^, an apparent horizon satisfies the equation 

H = = ViTi' + Kiju'n^ -K = 0, (1) 

where ra* is the outward-pointing unit normal to the horizon, all the field variables are 
evaluated on the horizon surface, and where for future use we define the "horizon function" 
H = ^"^^H as the left hand side of Eq. (|ip. (Notice that in order for the 3-divergence V^n* 
to be meaningful, ra* must be (smoothly) continued off the horizon, and extended to a field 
(3)^1 some 3-dimensional neighborhood of the horizon. The off- horizon continuation is 
non- unique, but it's easy to see that this doesn't affect H on the horizon.) 

To solve the apparent horizon equation Eq. (^, we begin by assuming that the horizon 
and coordinates are such that each radial coordinate line {(^,0) = constant} intersects the 
horizon in exactly one point. In other words, we assume that the horizon's coordinate shape 
is a "Strahlkorper" , defined by Minkowski as "a region in n-dimensional Euclidean space 
containing the origin and whose surface, as seen from the origin, exhibits only one point 
in any direction" (Ref. f^). Given this assumption, we can parameterize the horizon's 
shape by r = h{6, (p) for some single-valued "horizon shape function" h defined on the 
2-dimensional domain A of angular coordinates [6, 0). 

Equivalently, we may write the horizon's shape as ^^^F = 0, where the scalar function 
defined on some 3-dimensional neighborhood J\f of the horizon, satisfies ^^'^F = if and 
only if r = h{9, 0), and we take ^^'^F to increase outward. In practice we take ^^^F{r, 9, 0) = 
r-h{e,(f)). 

We define the non-unit outward-pointing normal (field) to the horizon by 

5. = = V.(%, (2) 

i.e. by 

Sr = 1 (3a) 
Su = -duh . (3b) 

and the outward-pointing unit normal (field) to the horizon by 



5 



Henceforth we drop the prefixes on ^^^Sj and ^^^n*. 

Substituting Eq. (H) into the apparent horizon equation Eq. (0), we see that the horizon 
function H{h) depends on the (angular) 2nd derivatives of h. In fact, the apparent horizon 
equation Eq. (|l]) is a 2nd order elhptic (boundary-value) PDE for h on the domain of angular 
coordinates A. The apparent horizon equation Eq. (|I|) must therefore be augmented with 
suitable boundary conditions to define a (locally) unique solution. These are easily obtained 
by requiring the horizon's 3-dimensional shape to be smooth across the artificial boundaries 
^ = 0, ^ = TT, = 0, and = 27r. 



IV. ALGORITHMS FOR SOLVING THE APPARENT HORIZON EQUATION 

We now survey various algorithms for solving the apparent horizon equation Eq. (P. 
Ref. reviews much of the previous work on this topic. 

In spherical symmetry, the apparent horizon equation Eq. (|I|) degenerates into a 1- 
dimensional nonlinear algebraic equation for the horizon radius h. This is easily solved 
by zero-finding on the horizon function H{h). This technique has been used by a number 
of authors, for example Refs. PD| , P^ , |52| ,|5[] . (See also Ref. |]r3| for an interesting analytical 



study giving necessary and sufficient conditions for apparent horizons to form in non- vacuum 
spherically symmetric spacetimes.) 

In an axisymmetric spacetime, the angular-coordinate space A is effectively 1-dimen- 
sional, so the apparent horizon equation Eq. (|ID reduces to a nonlinear 2-point boundary 
value (ODE) problem for the function h{6), which may be solved either with a shooting 
method, or with one of the more general methods described below. Shooting methods have 
been used by a number of authors, for example Refs. |r7| , |25| ,P], PT] , P!^ , |^ ,|^ . 



The remaining apparent-horizon-finding algorithms we discuss are all applicable to either 
axisymmetric spacetimes (2-dimensional codes) or fully general spacetimes (3-dimensional 
codes). 

Tod (Ref. [^]) has proposed an interesting pair of "curvature fiow" methods for finding 
apparent horizons. Bernstein (Ref. |jlO[) has tested these methods in several axisymmetric 
spacetimes, and reports favorable results. Unfortunately, the theoretical justification for 
these methods' convergence is only valid in time-symmetric {Kij = 0) slices. 

The next two algorithms we discuss are both based on a pseudospectral expansion of 
the horizon shape function h{6, 0) in some complete set of basis functions (typically spher- 
ical harmonics or symmetric trace- free tensors), using some finite number of the expansion 
coefficients {a^} to parameterize of the horizon shape. One algorithm rewrites the appar- 
ent horizon equation H{ak) = as ||iJ(afc)|| = 0, then uses a general-purpose function- 
minimization routine to search {afcj-space for a minimum of \\H\\. This algorithm has been 



used by Refs. |T5|;^ in axisymmetric spacetimes, and more recently by Ref. |^ in fully 
general spacetimes. Alternatively, Nakamura, Oohara, and Kojima (Refs. |^,48,47]) have 



suggested a functional iteration method for directly solving the apparent horizon equation 
H{ak) = for the expansion coefficients {a^}, and have used it in a number of fully general 
spacetimes. Kemball and Bishop (Ref. [0) have suggested and tested several modifications 
of this latter algorithm to improve its convergence properties. 

The final algorithm we discuss, and the main subject of this paper, poses the apparent 
horizon equation H{h) = as a nonlinear elliptic (boundary-value) PDE for h on the 
angular-coordinate space A. Finite differencing this PDE on an angular-coordinate grid 
{(^K; 0k)} gives a set of simultaneous nonlinear algebraic equations for the unknown values 
{/i(^k,0k)}, which are then solved by some variant of Newton's method. This "Newton's- 
method" algorithm (we continue to use this term even if a modification of Newton's method is 
actually used) has been used in axisymmetric spacetimes by a number of authors, for example 
Refs. [^,^,^,^,0 , and is also applicable in fully general spacetimes when the coordinates 
have a (locally) polar spherical topology. Huq (Ref. [0) has extended this algorithm to fully 
general spacetimes with Cartesian-topology coordinates and finite differencing, and much of 
our discussion remains applicable to his extension. 

The Newton's-method algorithm has three main parts: the computation of the discrete 
horizon function H(h), the computation of the discrete horizon function's Jacobian matrix 
J H(h) , and the solution of the simultaneous nonlinear algebraic equations H(h) = 0. We 
now discuss these in more detail. 



V. COMPUTING THE HORIZON FUNCTION 

In this section we discuss the details of the computation of the discrete horizon func- 
tion H(h). More precisely, first fix an angular-coordinate grid {(^K;0k)}- Then, given a 
"trial horizon surface" r = h{6,(j)), which need not actually be an apparent horizon, we 
define h{9, 0) to be the discretization of h{9, 0) to the angular-coordinate grid, and we con- 
sider the computation of H(h) on the discretized trial horizon surface, i.e. at the points 

{(r = h(^K,0K),^ = ^K,0 = 0K)}. 

The apparent horizon equation Eq. (|l]) defines H = ^"^^H in terms of the field variables 
and their spatial derivatives on the trial horizon surface. However, these are typically known 
only at the (3-dimensional) grid points of the underlying 3 + 1 code of which the horizon 
finder is a part. We therefore extend ^"^^H to some (3-dimensional) neighborhood J\f of the 
trial horizon surface, i.e. we define an extended horizon function ^^^H on J\f, 

(3)if = Viri^ + KijTi'n^ - K (7) 
= diu' + {di In ^)n' + Ki^n^n^ -K. (8) 

To compute *^^''H(h) on the (discretized) trial horizon surface, we first compute '^^)H(h) on 
the underlying 3 + 1 code's (3-dimensional) grid points in A/", then radially interpolate these 
^^^H values to the trial-horizon-surface position to obtain (^''H(h), 

(2)h(^,0) =interp((3)H(r,^,0),r = h(^,0)) , (9a) 

or equivalently 
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(2) 



Hi = interp((3)H<,i),r = hi 



(9b) 



where i is an angular grid-point index and the (ri) subscript denotes that the interpo- 
lation is done independently at each angular coordinate along the radial coordinate line 
{6 = 6i,(f) = (pi}. In practice any reasonable interpolation method should work well here: 
Refs. m. 



report satisfactory results using a spline interpolant; in this work we use a 
Lagrange (polynomial) interpolant centered on the trial-horizon-surface position, also with 
satisfactory results. Neglecting the interpolation error, we can also write Eq. (H) in the form 



{2)H( 



(3)H( 



9, 



(10) 



We consider two basic types of methods for computing the extended horizon function 
(3)H(h): 

• A "2-stage" computation method uses two sequential numerical finite differencing 
stages, first explicitly computing Sj and/or n* by numerically finite differencing h, 
then computing ^^^^H by numerically finite differencing Sj or n*. 

• A "1-stage" computation method uses only a single numerical (2nd) finite differencing 
stage, computing '^^^H directly in terms of h's 1st and 2nd angular derivatives. 

Figure |I| illustrates this. 

To derive the detailed equations for these methods, we substitute Eqs. (^ and Eqs. (^ 
into Eq. (|): 



= diU^ + {di In 



9,; 



A B 

+ 



K 



+ 



-K 



C 



-K . 



2^3/2 ^1/2 ^ 

where the subexpressions A, B, C, and D are given by 



C 

D = g'^SiSj , 



(11) 
(12) 

(13) 
(14) 



(15a) 

(15b) 
(15c) 
(15d) 



I.e. 



A = {g""' - g'^'^d^hW^ - g^^d^h)du.h 



B 
C 



- W - g'^'d^h) [d,g'^' - 2{d,gnduh + {d,gn{duh){d,h) 
d..g" - (5.^7*")5„/i] - + {d, In V^) ((?''■ - g'^dji) 



D = g" - Ig'^dji + g^\dji)[d,h) . 



(16a) 

(16b) 
(16c) 
(16d) 
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Comparing the 1-stage and 2-stage methods, the 2-stage methods' equations are some- 
what simpler, so these methods are somewhat easier to implement and somewhat cheaper 
(faster) to compute. However, for a proper comparison the cost of computing the horizon 
function must be considered in conjunction with the cost of computing the horizon function's 
Jacobian. Compared to the 1-stage method, the 2-stage methods double the effective radius 
of the net H(h) finite differencing molecules, and thus have 2(4) times as many nonzero 
off-diagonal Jacobian elements for a 2(3)-dimensional code. In practice the cost of comput- 
ing these extra Jacobian elements for the 2-stage methods more than outweighs the slight 
cost saving in evaluating the horizon function. We discuss the relative costs of the different 
methods further in section IVID. 



VI. COMPUTING THE JACOBIAN 

In this section we discuss the details of the computation of the Jacobian matrix J H(h) 
of the horizon function H(h) on a given trial horizon surface. 



A. Computing the Jacobian of a Generic Function P(Q) 



We consider first the case of a generic function P{Q) in d dimensions, finite differenced 
using A^-point molecules. We define the Jacobian matrix of the discrete P(Q) function by 



P(Q) 



or equivalently by the requirement that 



_ dPi 



5Pt = 



P(Q + 5Q)-P(Q) 



P(Q) . -^Qj 



ij 



(17a) 



(17b) 



for any infinitesimal perturbation 5Q of Q. 

We assume that P is actually a local grid function of Q, so the Jacobian matrix is sparse. 
(For example, this would preclude the nonlocal 4th order "compact differencing" methods 
described by Refs. pT|,|36|.) We assume that by exploiting the locality of the discrete P(Q) 



function, any single Pj can be computed in 0(1) time, independent of the grid size. 



1. Computing Jacobians by Numerical Perturbation 



We consider two general methods for computing the Jacobian matrix J P(Q) . The first 
of these is the "numerical perturbation" method. This involves numerically perturbing Q 
and examining the resulting perturbation in P(Q), 



P(Q) 



IJ 



P(Q + yue 



P(Q) 



(18) 



where e^'^-' is a Kronecker-delta vector defined by 
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.=1 







1 if I = J 
otherwise 



(19) 



and /i is a "small" perturbation amplitude. This computation of the Jacobian proceeds by 
columns: for each j, Qj is perturbed, and the resulting perturbation in P(Q) gives the Jth 
column of the Jacobian matrix. 

The perturbation amplitude fi should be chosen to balance the truncation error of the 
one-sided finite difference approximation Eq. ( plSf ) against the numerical loss of significance 



caused by subtracting the nearly equal quantities P(Q + /ie'^'^^) and P(Q). Refs. P5|J53| 
discuss the choice of fi, and conclude that if P(Q) can be evaluated with an accuracy of 
e, then fi ^ ^/e "seems to work the best". In practice the choice of fi isn't very critical 
for horizon finding. Values of 10~^ to 10~^ seem to work well, and the inaccuracies in the 
Jacobian matrix resulting from these values of /i don't seem to be a significant problem. 

This method of computing Jacobians requires no knowledge of the P(Q) function's inter- 
nal structure. In particular, the P(Q) function may involve arbitrary nonlinear computations, 
including multiple sequential stages of finite differencing and/or interpolation. This method 
is thus directly applicable to the ^^)H(h) computation. 

Assuming that P(Q) is already known, computing J P(Q) by numerical perturbation 

requires a total of A^*^ Pi evaluations at each grid point, i.e. it requires a perturbed-Pi 
evaluation for each nonzero Jacobian element. 



2. Computing Jacobians by Symbolic Differentiation 



An alternate method of computing the Jacobian matrix J P(Q) is by "symbolic differ- 
entiation" . This method makes explicit use of the finite differencing scheme used to compute 
the discrete P(Q) function. 

Suppose first that the continuum P{Q) function is a position-dependent local linear 
differential operator, discretely approximated by a position-dependent local finite difference 
molecule M, 



Pi= ^ M(i)mQi+m 

MgM(I) 



(20) 



Differentiating this, we have 



P(Q) 



^dP^_i M(i)j_i if j-i e M(i) 
iJ ~ 9Qj I otherwise 



(21) 



so that the molecule coefficients at each grid point give the corresponding row of the Jacobian 
matrix. 

More generally, suppose P is a position-dependent local nonlinear algebraic function of 
Q and some finite number of Q's derivatives, say 



P = PiQ,d,Q,d,jQ). 



(22) 



Logically, the Jacobian matrix J P(Q) is defined (by Eq. (}17\ )) in terms of the lineariza- 
tion of the discrete (finite differenced) P(Q) function. However, as illustrated in figure |^, if 
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the discretization (the finite differencing scheme) commutes with the hnearization, we can 
instead compute the Jacobian by first hnearizing the continuum P{Q) function, then finite 
differencing this (continuum) hnearized function. (This method of computing the Jaco- 
bian is essentially just the "Jacobian part" of the Newton-Kantorovich algorithm for solving 
nonlinear elliptic PDEs.) 

That is, we first linearize the continuum P{Q) function. 



6P 



dQ ^^d{d,Q) 



6diQ + 



dP 



9P 



dP 



d{diQ) 



d,6Q 



dP 



5dijQ 



dij6Q . 



(23) 
(24) 



We then view the linearized function 6P{6Q) as a linear differential operator, and discretely 
approximate it by the position-dependent finite difference molecule 



- dP 



(25) 



where I is the identity molecule and dj and dij are finite difference molecules discretely 
approximating di and dij respectively. Finally, we apply Eq. (|2T|) to the molecule M defined 



by Eq. (^5]) to obtain the desired Jacobian matrix J 
In practice, there's no need to explicitly form t 



P(Q)J- 

le molecule M 



the Jacobian matrix 



elements can easily be assembled directly from the known I, dj, and djj molecule coefficients 
and the "Jacobian coefficients" dP/dQ, dP/d{diQ), and dP/d{dijQ). Once these coefficients 
are known, the assembly of the actual Jacobian matrix elements is very cheap, requiring only 
a few arithmetic operations per matrix element to evaluate Eqs. ( p5D and (pT]). The main 
cost of computing a Jacobian matrix by symbolic differentiation is thus the computation 
of the Jacobian coefficients themselves. Depending on the functional form of the P{Q) 
function, there may be anywhere from 1 to 10 coefficients, although in practice these often 
have many common subexpressions. 

In other words, where the numerical perturbation method requires a Pi evaluation per 
nonzero Jacobian element, the symbolic differentiation method requires the computation of 
"a few" Jacobian-coefficient subexpressions per Jacobian row. More precisely, suppose the 
computation of all the Jacobian coefficients at a single grid point is J times as costly as a 
Pi evaluation. Then the symbolic differentiation method is approximately N'^/ J times more 
efficient than the numerical perturbation method. 



B. Semantics of the Horizon Function Jacobian 



We now consider the detai 
Jacobian of H(h) = (2)H(h), J 



ed semantics of the horizon function Jacobian. We define the 



H(h) 
J 



= J 



(2)H(h) 



(2)H(h)J, by 



iJ dhj ' 

(where i and j are angular grid-point indices), or equivalently by the requirement that 



(26a) 
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5(2)Hi= [(2)H(h + 5h)-(2)H(h)l = j[(2)H(h) 



ij 



Sh, 



(26b) 



for any infinitesimal perturbation 5h. Here i and J are both angular (2- dimensional) grid- 
point indices. Notice that this definition uses the total derivative d^'^^H/dh. This is because 
*^^)H(h) is defined to always be evaluated at the position r = h(6',0) of the trial horizon 



surface, so the Jacobian J (^)H(h) 
at a fixed position due to a pertur 



must take into account not only the direct change in ^"^^ H 
nation in h, but also the implicit change in '^^-'H caused by 



the field-variable coefficients in '^^^H being evaluated at a perturbed position r = h(6',0). 



It's also useful to consider the Jacobian J 
^^)H(h), which we define analogously by 



(3)H(h) 



(3)H(h) 



IJ 



~dh7' 



of the extended horizon function 



(27a) 



or equivalently by the requirement that 

^{3)Hj= [(3)H(h + 5h) -(=^)H(h) 



(3)H(h) 



IJ 



(27b) 



for any infinitesimal perturbation 5h. Here i is a 3-dimensional grid-point index for ^^-'H, 
while J is an (angular) 2-dimensional grid-point index for h. In contrast to J *^^)H(h) , 

this definition uses the partial derivative d^^^H/dh. This is because we take (3''H(h) to be 
evaluated at a fixed position (a grid point in the neighborhood JV of the trial horizon surface) 
which doesn't change with perturbations in h, so J ^3)H(h) need only take into account the 
direct change in '^'^^H at a fixed position due to a perturbation in h. 



'•^-'H(h) thus has much simpler semantics than J *^^''H(h) 



We have found J 



very useful, both as an intermediate variable in the computation of J 
the next section), and also conceptually, as an aid to thinking about the Jacobians. 



(3)H(h) 
^^^H(h) (described in 



C. Computing the Horizon Function Jacobian 



Table H (discussed further in section |VI D|) summarizes all the Jacobian-computation 
methods in this paper, which we now describe in detail. We tag each method with a short- 
hand "code", which gives the method's basic properties: whether it computes J (^^H(h) 

directly or computes J (3^H(h) as an intermediate step, whether it uses symbolic differenti- 
ation or numerical perturbation, and whether it uses a 1-stage or a 2-stage horizon function 
computation. 

The simplest methods for computing J *^^)H(h) 

directly with (2)H(h) in angular-coordinate space, without computing j K'^^H(h) 

mediate step. Since *^^)H(h) isn't given by a simple molecule operation of the form Eq. (pOD, 
symbolic differentiation isn't directly applicable here. However, numerical perturbation in 
angular-coordinate space is applicable, using either a 1-stage or a 2-stage method to com- 
pute (^^H(h). We refer to the resulting Jacobian computation methods as the "2d.np.ls" 
and "2d.np.2s" methods respectively. 



are the "2-dimensional" ones, which work 

as an inter- 
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Our remaining methods for computing J '^^)H(h) are all "3-dimensional" ones, which first 

explicitly compute J '^^)H(h) , then compute J '^^)H(h) from this in the manner described 
below. 

If '^^^H(h) is computed using the 1-stage method, i.e. via Eqs. (|T^) and ([TBD , then ei- 
ther numerical perturbation or symbolic differentiation may be used to compute J *^^)H(h) . 
We refer to these as the "3d.np.ls" and "3d.sd.ls" methods respectively. The symbohc- 
differentiation Jacobian coefficients for the 3d.sd.ls method are tabulated in appendix 

Alternatively, if '^^)H(h) is computed using a 2-stage method, then J (^)H(h) may be 

computed either by the simple numerical perturbation of '^^)H(h) (the "3d.np.2s" method), or 
by separately computing the Jacobians of the individual stages and matrix-multiplying them 
together. For the latter case, either numerical perturbation or symbolic differentiation may 
be used to compute the individual-stage Jacobians, giving the "3d.np2.2s" and "3d.sd2.2s" 
methods respectively. The symbolic-differentiation Jacobian coefficients for the 3d.sd2.2s 
method are tabulated in appendix 0. 

For any of the 3-dimensional methods, once J *^^-*H(h) is known, we compute J '^^•'H(h) 
as follows: 



(2)H(h) 



ij 



71, VI 



rfh(^^J,0j) 

rf(3)H(r = h(^i,, 



),^i, 



rfh(^j,0j 
9(3)H(r, ^1,01 



ah(^j 



+ 



(by Eq. (10)) 
9(3)H(r, ^1,0] 



interp(j[(3)H(h) 



r=h(ei,0i) 

r = 



(rl>j' 



hi) +interp'((3)H(rf),r = hj 



(28) 
(29) 
(30) 

(31) 
(32) 



where the (ri) subscripts in Eq. (|32D denote that the interpolations are done along the radial 
line {9 = 9i,(f) = (pi}, analogously to Eq. (|^), and where we neglect the interpolation errors 
in Eq. (H). 

Notice that the interp'(. . .) term in Eq. (p^ ) may be computed very cheaply using 
the same '•^•'H data values used in computing ^^^H, cf. Eq. (^. (The number of '•^^H data 
points used in the radial interpolation at each angular grid position will probably have to be 
increased by one to retain the same order of accuracy in the interp'(. 



as in the interp(. . .) term.) It's thus easy to compute J *^^^H(h) 



once J 



term in Eq. (32) 
(3)H(h) 



is known. 



D. Comparing the Methods 

Table | summarizes all the horizon-function and Jacobian computation methods de- 
scribed in sections and VI C . The table also shows which Jacobian matrices the meth- 



ods use, the methods' measured relative CPU times in our axisymmetric-spacetime (2- 
dimensional) code (discussed further in appendix 0), and our estimates of the methods' 
approximate implementation effort (programming complexity). 
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As can be seen from the table, for our implementation the 3d.sd.ls method is by far the 
most efficient of the Jacobian computation methods, being about a factor of 5 faster than any 
of the numerical perturbation methods. In fact, the computation of the Jacobian J (^)H(h) 
by the 3d.sd.ls method is only 1.5-2 times more expensive than the simple evaluation of the 
horizon function '•^^H(h). 

The relative performance of the different methods will of course vary considerably 
from one implementation to another, and especially between axisymmetric-spacetime (2- 
dimensional) and fully-general-spacetime (3-dimensional) codes. However, counting the 
number of operations needed for each method shows that the 3d.sd.ls method should re- 
main the fastest for any reasonable implementation. (We omit details of the counting in 
view of their length and lack of general interest.) Notably, the 3d.sd.ls method's relative 
advantage over the other methods should be approximately a factor of the molecule diame- 
ter larger for fully-general-spacetime (3-dimensional) codes than for axisymmetric-spacetime 
(2-dimensional) codes such as ours. 

Considering now the implementation efforts required by the various methods, in general 
we find that these depend more on which Jacobian matrices are involved than on how 
the Jacobians are computed: The 2-dimensional methods, involving only J 

the easiest to implement, while the 3-dimensional methods involving (only) J 



(2)H(h) 
(2)H(h) 



are 



and 



*^^)H(h) are somewhat harder to implement. The 3-dimensional methods involving the 



individual-stage Jacobians J Sj(h) , J n*(h) 



(3)H( 



, and/or J 



(3)H( 



are considerably 



more difficult to implement, due to these Jacobians' more complicated sparsity patterns. 

All the Jacobian matrices are highly sparse, and for reasonable efficiency it's essential to 
exploit this sparsity in their storage and computation. We have done this in our code, and 
our CPU-time measurements and implementation-effort estimates all reflect this. We briefly 
describe our sparse- Jacobian storage scheme in appendix 0. This scheme is very efficient, 
but its programming is a significant fraction of the overall Jacobian implementation effort, 
especially for the individual- stage Jacobians. 

Comparing numerical perturbation and symbolic differentiation methods, we had pre- 
viously suggested (Ref. |^) that symbolic-differentiation Jacobian computations would be 
very difficult to implement, necessarily requiring substantial support from a (computer) 
symbolic computation system. Several colleagues have expressed similar opinions to us. We 
had also previous suggested (Ref. [^) that due to the structure of the H{h) function, a 
Jacobian-coefficient formalism of the type described in sections VI A 2 and |V1 C| would not be 
valid for the horizon function, so symbolic differentiation methods would require explicitly 
differentiating the finite difference equations. 

These suggestions have proven to be incorrect: using the Jacobian-coefficient formalism 
described in sections |VI A2| and |VI C| , only the continuum equations need be differentiated, 
and this is easily done by hand. More generally, using this formalism we find the actual pro- 
gramming of the symbolic differentiation methods to be only moderately more difficult than 
that of the numerical perturbation methods. Some of the Jacobian coefficients tabulated in 
appendix ^ are fairly complicated, but no more so than many other computations in 3 + 1 
numerical relativity. 

In order to be confident of the correctness of any of the Jacobian-computation methods 
except the simple 2-dimensional numerical perturbation ones, we feel that it's highly desir- 
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able to program an independent method (which may be programmed for simphcity at the 
expense of efficiency) and make an end-to-end comparison of the resulting Jacobian matri- 
ces. (We have successfully done this for each of the Jacobian matrices computed by each of 
the methods listed in table |, and our implementation-effort estimates there include doing 
this.) If, and only if, the Jacobians agree to within the expected truncation error of the 
numerical-perturbation Jacobian approximation Eq. (|l^), then we can have a high degree 
of confidence that both calculations are correct. If they disagree, then we find the detailed 
pattern of which matrix elements differ to be a very useful debugging aid. 

Summarizing our comparisons, then, we find that the best Jacobian computation method 
is clearly the 3d.sd.ls one. It's much more efficient than any of the other methods, and still 
quite easy to implement. 



VII. CONVERGENCE TESTS 

Before continuing our discussion of Newton's-method horizon finding, in this section we 
digress to consider the convergence of finite differencing computations to the continuum 
limit. 

As has been forcefully emphasized by Choptuik (Refs. |lT8|-pO[]), a careful comparison 
of a finite differencing code's numerical results at different grid resolutions can yield very 
stringent tests of the code's numerical performance and correctness. In particular, such 
"convergence tests" can yield reliable numerical estimates of a code's external errors, i.e. of 
the deviation of the code's results from those that would be obtained by exactly solving the 
continuum equations. With, and only with, such estimates available, we can safely draw 
inferences about solutions of the continuum equations from the code's (finite-resolution) 
numerical results. 

To apply this technique in the horizon- finding context, suppose first that the (a) true 
(continuum) apparent horizon position h* is known. For a convergence test in this case, we 
run the horizon finder twice, using a 1:2 ratio of grid resolutions. As discussed in detail by 



Ref. [19|, if the code's numerical errors are dominated by truncation errors from nth order 



finite differencing, the numerically computed horizon positions h must satisfy 

h[Ax] = h* + (Ax)"/ + 0((Ax)"+2) (33a) 
h[Ax/2] = h* + (Aa;/2)"/ + 0((A,t)"+2) (33b) 

at each grid point, where h[Ax] denotes the numerically computed horizon position using 
grid resolution Ax, and / is an 0(1) smooth function depending on various high order 
derivatives of h* and the field variables, but not on the grid resolution. (We're assuming 
centered finite differencing here in writing the higher order terms as 0((Ax)"'"'"^), otherwise 
they would only be 0{{Ax)"'~^^).) Neglecting the higher order terms, i.e. in the limit of small 
Ax, we can eliminate / to obtain a direct relationship between the code's errors at the two 
resolutions, 

h[Ax/2]-h* ^ ^ 
h[Ax]-h* 2-' ^ 

which must hold at each grid point common to the two grids. 
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To test how well any particular set of (finite-resolution) numerical results satisfies this 
convergence criterion, we plot a scatterplot of the high- resolution errors h[Ax/2] — /i* against 
the low-resolution errors h [Ax] — h* at the grid points common to the two grids. If, and given 
the arguments of Ref. |jT9|, in practice only if, the error expansions Eq. (|33|) are valid with the 
higher order error terms negligible, i.e. if and only if the errors are indeed dominated by the 
expected rath order finite difference truncation errors, then all the points in the scatterplot 
will fall on a line through the origin with slope 1/2". 

Now suppose the true (continuum) apparent horizon position h* is unknown. For a 
convergence test in this case, we run the horizon finder three times, using a 1:2:4 ratio of 
grid resolutions. Analogously to the 2-grid case, we now have 

h[Ax] =h* + (Ax)"/ + 0((Ax)"+2) (35a) 
h[Ax/2] =h* + (Ax/2)"/ + 0((Ax)"+2) (35b) 
h[Ax/4] =h* + (Ax/4)"/ + 0((Ax)"+2) , (35c) 

at each grid point, with / again independent of the grid resolution. Again neglecting the 
higher order terms, we can eliminate both / and h* to obtain the "3-grid" convergence 
criterion 

h [Ax/2] -h [Ax/4] ^ l_ 
h[Ax] - h[Ax/2] 2" ^ ' 

which must hold at each grid point common to the three grids. We test this criterion using 
a scatterplot technique analogous to that for the 2-grid criterion Eq. (|3^. 

We emphasize that for a 3-grid convergence test of this type, the true continuum solution 
h* need not be known. In fact, nothing in the derivation actually requires h* to be the true 
continuum horizon position - it need only be the true continuum solution to some continuum 
equation such that the truncation error formulas Eq. (^) hold. We make use of this latter 
case in sections |V111 (J| and [IX B| to apply 3-grid convergence tests to intermediate Newton 



iterates (trial horizon surfaces) of our horizon finder. 

For both the 2-grid and the 3-grid convergence test, we find that the pointwise nature 
of the scatterplot comparison makes it significantly more useful than a simple comparison 
of gridwise norms. In particular, the scatterplot comparison clearly shows convergence 
problems which may occur only in a small subset of the grid points (for example near a 
boundary), which would be "washed out" in a comparison of gridwise norms. 

Notice also that the parameter n, the order of the convergence, is (should be) known 
in advance from the form of the finite differencing scheme. Thus the slope-1/2" line with 
which the scatterplot points are compared isn't fitted to the data points, but is rather an 
a priori prediction with no adjustable parameters. Convergence tests of this type are thus 
a very strong test of the validity of the finite differencing scheme and the error expansions 
Eq. (H) or Eq. (|3|). 

VIII. SOLVING THE NONLINEAR ALGEBRAIC EQUATIONS 

Returning to our specific discussion of horizon finding, we now discuss the details of 
using Newton's method or a variant to solve the simultaneous nonlinear algebraic equations 
H(h) = 0. 
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A. Newton's Method 



The basic Newton's-method algorithm is well known: At each iteration, we first linearize 
the discrete H(h) function about the current approximate solution h^^^ 



(37) 



denotes the Jacobian 



H(h('=) + 5h) = H(h(*^)) + j[H(hW)] ■ 5h + 0{\\Shf) 

where Sh now denotes a finite perturbation in h, and where J H(h*^'^)' 
matrix J H(h) evaluated at the point h = h*^^). We then neglect the higher order (nonlinear) 

terms and solve for the perturbation 5h^''^ such that H(h^'^^ + ^h^'^^) = 0. This gives the 
simultaneous linear algebraic equations 



H(hW)] -^hW = -H(h('=)) 
to be solved for dh^'^K Finally, we update the approximate solution via 



(38) 



(39) 



and repeat the iteration until some convergence criterion is satisfied. 

Notice that here we're using the word "convergence" in a very different sense from that of 
section |V11| here it refers to the "iteration-convergence" of the Newton iterates h*^^-* to the 
exact solution h* of the discrete equations, whereas there it refers to the "finite-difference- 
convergence" of a finite difference computation result h [Ax] to its continuum limit h* as the 
grid resolution is increased. 

Once the current solution estimate h'^'^^ is reasonably close to h*, i.e. in practice once 
the trial horizon surface is reasonably close to the (an) apparent horizon, Newton's method 
converges extremely rapidly. In particular, once the linear approximation in Eq. ( ^7\) is 
approximately valid, Newton's method roughly squares the relative error ||h — h*||/||h*|| 
at each iteration, and can thus bring the error down to a negligible value in only a few 
(more) iterations. (This rapid "quadratic" convergence depends critically on the mutual 
consistency of the horizon function and Jacobian matrix used in the computation, and is 
thus a useful diagnostic for monitoring the Jacobian's correctness.) (For a detailed discussion 
of Newton's method, including precise formulations and proofs of these statements, see, for 
example, Ref. pl[.) 

However, if the initial guess h'^°^ for the horizon position, or more generally any Newton 
iterate (trial horizon surface) h*^'^^ differs sufficiently from h* so that the linear approximation 
in Eq. ( pTj ) isn't approximately valid, then Newton's method may converge poorly, or fail to 
converge at all. 



B. Modifications of Newton's Method 



Unfortunately, as discussed in section IX B, for certain types of initial guesses Newton's 



method fails to converge unless the initial guess is very close to the exact solution of the 
finite difference equations. There's an extensive numerical analysis literature on more ro- 
bust "modified Newton" algorithms for solving nonlinear algebraic equations, for example 
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Refs. [p|,p|, p6| , ^ , p0| , |66| , |67| . We have found Ref. p6| to be a particularly useful introduction 
to this topic. 

For horizon finding, the Jacobian matrix's size is the number of angular grid points on 
the horizon surface. This is generally large enough that it's important for the nonlinear- 
algebraic-equations solver to support treating the Jacobian as either a band matrix (for 
axisymmetric-spacetime codes) or a fully general sparse matrix (for fully-general-spacetime 
codes). It's also desirable for the nonlinear- algebraic-equation solver to permit explicit 
bounds on the solution vector, so as to ensure the trial horizon surfaces never fall outside 
the radial extent of the code's main 3-dimensional grid. Unfortunately, these requirements 
rule out many nonlinear-algebraic-equation software packages. 

For the sake of expediency, in the present work we chose to write our own implemen- 
tation of a relatively simple modified-Newton algorithm, the "line search" algorithm de- 
scribed by Refs. p6| , |50| . However, a much better long-term solution would be to use an 
extant nonlinear-algebraic-equations code embodying high-quality implementations of more 
sophisticated algorithms, such as the GIANT code described by Refs. [^^. We would 
expect Newton's-method horizon-finding codes using such software to be considerably more 
robust and efficient than our present code. 

The modified-Newton algorithm used in this work, the line-search algorithm of 
Refs. is identical to the basic Newton's-method algorithm, except that the Newton's- 

method update Eq. ([39|) is modified to h^'^^^) ^ h'^'^^ + XSh^''\ where A G (0,1] is chosen 
at each "outer" iteration by an inner "line search" iteration to ensure that ||H||2 decreases 
monotonically. Refs. p6| , |50[| show that such a choice of A is always possible, and describe an 
efficient algorithm for it. Sufficiently close to the solution h*, this algorithm always chooses 
A = 1, and so takes the same steps as Newton's method. The overall modified-Newton 
algorithm thus retains the extremely rapid convergence of Newton's method once the linear 
approximation in Eq. (^) is good. 

The line-search algorithm described by Refs. [^,0] always begins by trying the basic 
Newton step A = 1. For horizon finding, we have slightly modified the algorithm to decrease 
the starting value of A if necessary to ensure that h^^^ + A Sh^'^^ lies within the radial extent 
of our code's main (3-dimensional) numerical grid at each angular grid coordinate. Our im- 
plementation of the algorithm also enforces an upper bound (typically 10%) on the relative 
change ||A 5h'-'^''/h'~'^''||oo in any component of h^^^ in a single outer iteration. However, it's not 
clear whether or not this latter restriction is a good idea: although it makes the algorithm 
more robust when the H(h) function is highly nonlinear, it may slow the algorithm's conver- 
gence when the H(h) function is only weakly nonlinear and the error in the initial guess is 
large. We give an example of this latter behavior in section ^ 



C. The Newton-Kantorovich Method 

We have described the Newton's-method algorithm, and the more robust modified ver- 
sions of it, in terms of solving the discrete H(h) = equations. However, these algorithms 
can also be interpreted directly in terms of solving the continuum H{h) = equations. This 
"Newton-Kantorovich" method, and its relationship to the discrete Newton's method, are 



discussed in detail by Ref. |14 . 



18 



For the Newton-Kantorovich algorithm, at each iteration, we first hnearize the continuum 
differential operator H{h) about the current continuum approximate solution h^^\ 



H{h^^^ + 5h) = H{h^^^) + {6h) + 0{\\6h\\ 



(40) 



where Sh is now a finite perturbation in h, and where the linear differential operator 
J H{h'^^^) is now the linearization of the differential operator H{h) about the point h = h'^^\ 

We then neglect the higher order (nonlinear) terms and solve for the perturbation 8h'^^^ such 
that H{h^^'^ + 5/i(^)) = 0. This gives the linear differential equation 



J 



to be solved for Sh^'^^ Finally, we update the approximate solution via 



(41) 



(42) 



and repeat the iteration until some convergence criterion is satisfied. 

Now suppose we discretely approximate this continuum Newton-Kantorovich algorithm 
by finite differencing the iteration equation Eq. (^). If the finite differencing and the lin- 
earization commute in the manner discussed in section [VI A 2| , then this finite-difference 
approximation to the Newton-Kantorovich algorithm is in fact identical to the discrete 
Newton' s-method algorithm applied to the (discrete) H(h) = equations obtained hy finite dif- 
ferencing the continuum H{h) = equation. (In a simpler context, our Jacobian-coefficient 
formalism described in section [VIA 2 essentially just exploits the "Jacobian part" of this 
identity.) 

Therefore, when using the discrete Newton's method to solve the H(h) = equations, we 
can equivalently view each Newton iterate (trial horizon surface) h*^^)[Aa;] as being a finite 
difference approximation to the corresponding continuum Newton-Kantorovich iterate (trial 
horizon surface) h^''\ As the grid resolution is increased, each Newton iterate h*^'^)[Aa;] should 
therefore show proper finite-difference-convergence regardless of the iteration- convergence or 
iteration- divergence of the Newton or Newton-Kantorovich iteration itself. 

Moreover, once we verify the individual Newton iterates' finite-differencing-convergence 
(with a 3-grid convergence test), we can safely extrapolate the iteration-convergence or 
iteration-divergence of this discrete iteration to that of the continuum Newton-Kantorovich 
algorithm applied to the (continuum) H{h) = equations. In other words, by this procedure 
we can ascribe the iteration-convergence or iteration-divergence of Newton's method to 
inherent properties of the continuum H{h) = equations, as opposed to (say) a finite 
differencing artifact. We make use of this in section |IXB. 



IX. GLOBAL CONVERGENCE OF THE HORIZON FINDER 

We now consider the global convergence behavior of the Newton's-method horizon finding 
algorithm. That is, how close must the initial guess h*^°^ be to the (an) exact solution h* of 
the finite difference equations in order for the iterates (trial horizon surfaces) h'^'^^ to converge 
to h*? In other words, how large is the algorithm's radius of convergence? 
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A. Global Convergence for Schwarzschild Spacetime 



To gain a general picture of the qualitative behavior of H{h) and its implications for 
Newton's-method horizon finding, it's useful to consider Schwarzschild spacetime. We use 
the Eddington-Finkelstein slicing, where the time coordinate is defined by requiring t + r 
to be an ingoing null coordinate. (These slices aren't maximal: K is nonzero and spatially 
variable throughout the slices.) 

Taking the black hole to be of dimensionless unit mass, the (only) apparent horizon in 
such a slice is the coordinate sphere r = 2. More generally, a straightforward calculation 
gives 

H = ^j'--JL (43) 

for spherical trial horizon surfaces with coordinate radius r. Figure |^ shows H{r) for these 
surfaces. As expected, H = for the horizon r = 2. However, notice that H reaches a 
maximum value at r = r™^^ = |(3 + VSS) ~ 4.372, and in particular that for r > r™*^^, 
H > and dH/dr < 0. Because of this, almost any algorithm - including Newton's method 
and its variants - which tries to solve H{r) = using only local information about H{r), 
and which maintains the spherical symmetry, will diverge towards infinity when started from 
within this region, or if any intermediate iterate (trial horizon surface) ever enters it. 

In fact, we expect broadly similar behavior for H in any black hole spacetime: Given 
an asymptotically flat slice containing an apparent horizon or horizons, consider any 1- 
parameter family of topologically 2-spherical nested trial horizon surfaces starting at the 
outermost apparent horizon and extending outward towards the 2-sphere at spatial infinity. 
H = for the horizon and for the 2-sphere at spatial infinity, so \\H\\ must attain a maximum 
for some finite trial horizon surface somewhere between these two surfaces. We thus expect 
the same general behavior as in the Schwarzschild-slice case, i.e. divergence to infinity if the 
initial guess or any intermediate iterate (trial horizon surface) lies outside the maximum- 
surface. This argument isn't completely rigorous, since the algorithm could move inward in 
an angularly anisotropic manner, but this seems unlikely. 

Fortunately, in practice this isn't a problem: the black hole area theorem places an upper 
bound on the size of an apparent horizon, and this lets us avoid overly-large initial guesses, 
or restart the Newton iteration if any intermediate iterate (trial horizon surface) is too large. 



B. Global Convergence in the Presence of High-Spatial-Frequency Errors 

Assuming the initial guess is close enough to the horizon for the divergence-to-infinity 
phenomenon not to occur, we find the global convergence behavior of Newton's method 
to depend critically on the angular spatial frequency spectrum of the initial guess's error 
h(o) _ h*; If j^iiQ error has only low-spatial-frequency components (in a sense to be clarified 
below), then Newton's method has a large radius of convergence, i.e. it will converge even for 
a rather inaccurate initial guess. However, if the error has significant high- spatial- frequency 
components, then we find that Newton's method has a very small radius of convergence, i.e. it 
often fails to converge even when the error 

h(o)_ 

h* is very small. 
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This behavior is not an artifact of insufficient resolution in the finite difference grid. 
Rather, it appears to be caused by a strong nonhnearity in the continuum H{h) function 
for high-spatial-frequency components in h. In this context there's no sharp demarcation 
between "low" and "high" spatial frequencies, but in practice we use the terms to refer to 
angular Fourier components varying as (say) cos m9 with m ^ 4 and m ^ 8 respectively. 

1. An Example 

As an example of this behavior, consider Kerr spacetime with dimensionless angular 
momentum a = J/M"^ = 0.6. We use the Kerr slicing, where the time coordinate is defined 
by requiring t + r to be an ingoing null coordinate. (These slices generalize the Eddington- 
Finkelstein slices of Schwarzschild spacetime, and are similarly nonmaximal, with K nonzero 
and spatially variable throughout the slices.) Taking the black hole to be of dimensionless 
unit mass, the (only) apparent horizon in such a slice is the coordinate sphere r = h*{6, cj)) = 
1 + = 1.8. 

For this example we consider two different initial guesses for the horizon position: one 
containing only low-spatial-frequency errors, r = h^o\9,(P) = 1.8 + 0.1 cos 4^, and one con- 
taining significant high-spatial-frequency errors, r = h^^\6,(f)) = 1.8 + 0.1 cos 10^. Notice 
that both initial guesses are quite close to the actual horizon shape, differing from it by 
slightly less than 5%. We use a finite difference grid with AO = which is ample to 
resolve all the trial horizon surfaces occurring in the example. 

Figure §(a) shows the behavior of Newton's method for the low-spatial- frequency-error 
initial guess. As can be seen, here Newton's method converges without difficulty. 

Figure ^(b) shows the behavior of Newton's method for the high-spatial- frequency-error 
initial guess. Here Newton's method fails to converge: the successive iterates (trial horizon 
surfaces) h^'^) move farther and farther away the horizon, and rapidly become more and more 
nonspherical. 

Figure ^(c) shows the behavior of the modified Newton's method for this same high- 
spatial-frequency-error initial guess. Although the first iteration still moves the trial horizon 
surface somewhat inward from the horizon, the nonsphericity damps rapidly, and the suc- 
cessive iterates (trial horizon surfaces) quickly converge to the horizon. 

Notice that all the intermediate iterates (trial horizon surfaces) in this example are well- 
resolved by the finite difference grid. To verify that insufficient grid resolution isn't a factor 
in the behavior of the horizon finder here, we have rerun all three parts of this example with 
several higher grid resolutions, obtaining results essentially identical to those plotted here. 

More quantitatively, following our discussion of the Newton-Kantorovich method in sec- 
tion [VIII Q we have made 3-grid convergence tests of each intermediate iterate (trial horizon 
surface) in this example. For example, figure ^ shows a 3-grid convergence test for the New- 
ton iterate (trial horizon surface) h'-^^ plotted in figure §(b), using grids with resolutions 
A9 = ^^Yof'iof ■ Notice that despite the iteration-divergence of the Newton iteration, 
this iterate shows excellent 4th order finite-difference-convergence. The other Newton and 
modified-Newton iterates (trial horizon surfaces) in our example all similarly show excellent 
4th order finite-difference-convergence. 

We conclude that the iteration-divergence of Newton's method seen in figure §(b), is in 
fact an inherent property of the continuum Newton-Kantorovich algorithm for this initial 
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guess and slice. Looking at the internal structure of this algorithm, we see that its only 
approximation is the linearization of the continuum H{h) function in Eq. (0), so the algo- 
rithm's iteration-divergence must (can only) be due to nonlinearity in the continuum H{h) 
function. 

2. The Horizon-Perturbation Survey 

To investigate how general the poor convergence of Newton's method seen in this example 
is, and to what extent it also occurs for the modified Newton's method, we have made a 
Monte Carlo numerical survey of both algorithms' behavior over a range of different initial- 
guess-error spatial frequency spectra. 

For this survey we first fix a particular horizon-finding algorithm. Suppose we are given 
a slice containing an apparent horizon at the continuum position h*, and consider running 
the horizon finder with the generic perturbed initial guess 

M 

h = h*+ Cm cos mO (44) 

m=0 
m even 

for some set of initial-guess-error Fourier coefficients {cm}- (Here we include only even-m 
cosine terms in 9 so as to preserve axisymmetry and equatorial reflection symmetry, which 
our code requires.) 

For each value of M we deflne the horizon flnder's "convergence region" in {cm}-space 
to be the set of coefficients {cm} for which the horizon finder converges (we presume to the 
correct solution). For example, the convergence region will in practice always include the 
origin in {cm}-space, since there h = h*, so the initial guess differs from the exact solution 
of the discrete H(h) = equations only by the small H(h) finite differencing error. 

We define Vm to be the (hyper)volume of the convergence region. As described in detail 
in appendix 0, we estimate Vm by Monte Carlo sampling in {cm,}-space. Given Vm, we then 
define the "volume ratio" 

(Vo if M = 
ifM>2 • 

[ Vm-2 

SO that Rm measures the average radius of convergence of the horizon finder in the cm 
dimension. 

3. Results and Discussion 

We have carried out such a horizon-perturbation survey for the same Kerr slices of the 
unit-mass spin-0.6 Kerr spacetime used in the previous example, for both the Newton and 
the modified-Newton algorithms, for M = 0, 2, 4, . . . , 12. Figure ^ shows the resulting 
volume ratios. Although the precise values are somewhat dependent on the details of our 
implementation and on the test setup (in particular on the position of the inner grid bound- 
ary, which is at r = 1 for these tests), the relative trends in the data should be fairly generic. 
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These tests use a grid with A6' = which is adequate to resolve all the perturbed trial 
horizon surfaces. 

As can be seen from the figure, the modified-Newton algorithm is clearly superior to the 
Newton algorithm, increasing the radius of convergence by a factor of 2-3 at high spatial 
frequencies. However, both algorithms' radia of convergence still fall rapidly with increasing 
spatial frequency, approximately as l/M^^^, although the rate is slightly slower for the 
modified-Newton than for the Newton algorithm. The radius of convergence of Newton's 
method falls below 0.1 (~ 5% of the horizon radius) by M ^ 10, and the data suggest 
that the radius of convergence of the modified-Newton method would be similarly small 
by M > 18. 

Since the grid resolution is adequate, we again conclude that the small radius of conver- 
gence of Newton's method must be due to a strong high-spatial-frequency nonlinearity in 
the continuum H{h) function. Our horizon-perturbation survey covers only a single axisym- 
metric initial slice and generic axisymmetric perturbations of the initial guess, but it seems 
unlikely that the nonlinearity would diminish for more general cases. Huq (Ref. |^) has 
made limited tests with nonaxisymmetric spacetimes and high-spatial-frequency perturba- 
tions, and has found (poor) convergence of Newton's method similar to our results. 

Although we write the continuum horizon function as H = H{h), it's more accurate to 
write this as H = H{gij, Kij, h), since H also depends on the slice's field variables and their 
spatial derivatives. Examining the functional form of the H{gij, K^j, h) function in Eqs. ([M]) 
and (|T^) , we see that H depends on the g^^ components in a manner broadly similar to its 
dependence on h. We thus conjecture that the H{gij, Kij, h) function may exhibit strong 
high-spatial-frequency nonlinearity in the field variables, in particular in the g^^ components, 
similar to its nonlinear dependence on h. 

If this is the case, then high-spatial-frequency variations in the field variables, such as 
would be caused by high-frequency gravitational radiation, might well impair the conver- 
gence of Newton's method in a manner similar to high-spatial-frequency perturbations in 
h. Further investigation of this possibility, either by analytical study of the nonlinear struc- 
ture of the H{gij, Kij, h) function, or by numerical investigations, would be very interesting. 
Fortunately, however, those (few) dynamic black hole spacetimes which have been explicitly 
computed thus far (for example Ref. 0) seem to contain mainly low-frequency gravitational 
radiation. 

In general, how serious a problem is the poor high-spatial-frequency convergence of New- 
ton's method? Given a sufficiently good initial guess, Newton's method still converges very 
rapidly (quadratically), so the key question is, how good is the initial guess in practice? Two 
cases seem to be of particular importance: If the horizon finder is being used to update a 
horizon's position at each time step of a time evolution, then the previous time step's horizon 
position probably provides a sufficiently good initial guess for Newton's method to converge 
well. In contrast, if the horizon finder is being used on initial data, or in a time evolution 
where there is no nearby horizon in the previous time step, then significant initial-guess 
errors can be expected, and Newton's method may converge poorly. 
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X. ACCURACY OF THE HORIZON FINDER 



We now consider the accuracy of the Newton's-method horizon finding algorithm. That 
is, assuming the Newton or modified-Newton iteration converges, how close is the hori- 
zon finder's final numerically computed horizon position to the (a) true continuum horizon 
position h*l 

The horizon finder computes Newton or modified-Newton iterates (trial horizon surfaces) 
|-|(A:) fQ]^ = 0, 1, 2, . . . , Until somc convergence criterion is satisfied, say at k = p. Because 
of the extremely rapid convergence of the Newton and modified-Newton iterations once the 
error is sufficiently small (cf. section |V111 A|), there's little extra cost in using a very strict 



convergence criterion, i.e. in solving the discrete H(h) = equations to very high accuracy. 
In our horizon finder we typically require || H(h^P)) ||oo < 10~^°. 

We denote the exact solution of the discrete H(h) = equations by h*. Given that 
||H(h*^P))|| is reasonably small, then from standard matrix-perturbation theory (see, for ex- 
ample, Refs. 123,|32|), ||h(p) - h*|| < K||H(h(p))||, where K is the condition number of the 
(presumably nonsingular) Jacobian matrix J H(h) at the horizon position. 

If we take the convergence tolerance to be strict enough for Hh*^^) — h*|| to be negligible, 
then the overall accuracy of the horizon finder, i.e. the external error ||h(p) - h*\\ in the 
computed horizon position, is thus limited only by the closeness with which the discrete 
H(h) =0 equations approximate the continuum H{h) = equations, i.e. by the accuracy 
of the H(h) finite differencing. This potential for very high accuracy is one of the main 
advantages of the Newton's-method horizon-finding algorithm. 

For an example of the accuracy attainable in practice, we again consider the Kerr slices 
of the unit-mass spin-0.6 Kerr spacetime. However, to make the horizon deviate from a 
coordinate sphere and hence be a more significant test case for our horizon finder, we apply 
the spatial coordinate transformation 

r r + — (a2 cos 26 + cos 46') (46a) 

5^ _l_ V y 

to the slice, where the parameters are given by 

b = 5 02 = 0.75 04 = 0.05. (46b) 

As shown in figure |3(a), in the transformed coordinates this gives a strongly non-spherical 
"peanut- shaped" horizon, similar in shape to those around a pair of coalescing black holes. 

We have run our horizon finder on this slice, using the warped-coordinate coordinate 
sphere r = 1.8 as an initial guess and a grid resolution of A9 = We used the modified- 
Newton algorithm, which converged to the horizon without difficulty. (The convergence took 
9 iterations, but would have taken only 6 iterations in the absence of our 10% restriction 
on the relative change in any component of h in a single outer iteration, cf. section |VlIIIj| .) 



Figure |^(a) shows the initial guess and the final numerically computed horizon position. 

Figure |^(b) shows the results of a 2-grid convergence test of the final numerically com- 
puted horizon position for this example, using grids with resolutions A9 = ^:Yoi' 
seen, the numerically computed solution shows excellent 4th order convergence. Moreover, 
the numerically computed horizon positions are very accurate, with \\h^P'^-h*\\ ~ 10-^(10"^) 
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for a grid resolution of = ^(yoI)- Errors of this magnitude are typical of what we find 
for Newton's-method horizon finding using 4th order finite differencing, so long as the grid 
adequately resolves the horizon shape. 



XI. FINDING OUTERMOST APPARENT HORIZONS 

The main focus of this paper is on locally finding apparent horizons, i.e. on finding an 
apparent horizon in a neighborhood of the initial guess. However, there's a related global 
problem of some interest which has heretofore attracted little attention, that of finding or 
recognizing the outermost apparent horizon in a slice. (By "recognizing" the outermost 
apparent horizon we mean the problem of determining whether or not a given apparent 
horizon is in fact the outermost one in a slice.) 

These global problems are of particular interest when apparent horizons are used to set 
the inner boundary of a black-hole-excluding grid in the numerical evolution of a multiple- 
black-hole spacetime, as discussed by Refs. ||57|j5^j58| j5[] . In this context, we can use the 



appearance of a new outermost apparent horizon surrounding the previously-outermost ap- 
parent horizons around two black holes as a diagnostic that the black holes have collided 
and coalesced into a single (distorted) black hole. As suggested by Ref. [^], we can then 



generate a new numerical grid and attach it to the new outermost apparent horizon, and 
continue the evolution on the exterior of the new (distored) black hole. 

So far as we know, no reliable algorithms are known for finding or recognizing outermost 
apparent horizons in nonspherical spacetimes. (For spherical spacetimes, a 1-dimensional 
search on H{r) suffices.) If started with a very large 2-sphere as the initial guess, the 
curvature flow method might well converge to the outermost horizon in the slice, but as 
mentioned in section the theoretical justification for this method's convergence is only 
valid in time-symmetric {Kij = 0) slices. 



For the remaining local-horizon-finding algorithms surveyed in section including the 
Newton's-method one, we know of no better method for locating or recognizing outermost 
horizons than trying the local-horizon-finder with a number of different initial guesses near 
the suspected position of an outermost horizon. If this method succeeds it locates a horizon, 
but there's still no assurance that this horizon is the outermost one in the slice. Moreover, if 
all the local-horizon-finding trials fail, this may mean that there's no horizon in the vicinity 
of the initial guesses, or it may only mean that a horizon is present nearby but the method 
failed to converge to it. It's also not clear how many local- horizon-finding trials should be 
made, nor just how their initial guesses should be chosen. 

This is clearly not a satisfactory algorithm. Further research to develop reliable al- 
gorithms for finding or recognizing outermost apparent horizons in generic (nonspherical, 
nonmaximal) slices would be very useful. 



XII. CONCLUSIONS 

We find Newton's method to be an excellent horizon-finding algorithm: it handles fully 
generic slices, it's fairly easy to implement, it's very efficient, it's generally robust in its 
convergence, and it's very accurate. These properties are all well known, and Newton's 
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method is widely used for horizon finding. In this paper we focus on two key aspects of this 
algorithm: the computation of the Jacobian matrix, and the algorithm's global convergence 
behavior. 

Traditionally, the Newton's-method Jacobian matrix is computed by a numerical pertur- 
bation technique. In this paper we present a much more efficient "symbolic differentiation" 
technique. Conceptually, this entails differentiating the actual finite difference equations 
used to compute the discrete horizon function H(h). However, provided the finite differ- 
encing scheme commutes with linearization, the computation can instead be done by first 
differentiating the continuum horizon function H{h), then finite differencing. (This is es- 
sentially just the "Jacobian part" of the Newton-Kantorovich method for solving nonlinear 
PDEs.) 

In our axisymmetric-spacetime (2-dimensional) numerical code, this method is about a 
factor of 5 faster than than any other Jacobian computation method. In fact, the Jacobian 
computation using this method is only 1.5-2 times more expensive than the simple evaluation 
of H(h). We expect the symbolic differentiation method's relative advantage over other 
Jacobian computation methods to be roughly similar for other axisymmetric-spacetime (2- 
dimensional) codes, and an additional factor of ~ 3-5 larger for fuUy-general-spacetime 
(3-dimensional) codes. 

We had previously suggested (Ref. |5^) that symbohc-differentiation Jacobian compu- 
tations would be quite difficult, necessarily requiring substantial support from a (computer) 
symbolic computation system. Several colleagues have expressed similar opinions to us. 
However, this turns out not to be the case: we computed all the symbolic-differentiation 
Jacobian coefficients for our horizon finder by hand in only a few few pages of algebra. Some 
of the coefficients are fairly complicated, but no more so than many other computations in 
3 + 1 numerical relativity. 

We find the actual programming of the symbolic differentiation Jacobian computation 
to be only moderately more difficult than that of a numerical perturbation computation. In 
order to be confident of the correctness of a symbolic differentiation Jacobian computation, 
we feel that it's highly desirable to program an independent numerical perturbation method 
and make an end-to-end comparison of the resulting Jacobian matrices. The comparison 
Jacobian computation may be programmed for simplicity at the expense of efficiency, so it 
needn't add much to the overall symbolic-differentiation implementation effort. 

Turning now to the convergence behavior of Newton's method, we find that so long as 
the error in the initial guess (its deviation from the true horizon position) contains only low- 
spatial-frequency components, a Newton's-method horizon finder has a large (good) radius 
of convergence, i.e. it converges even for rather inaccurate initial guesses. However, if the 
error in the initial guess contains significant high-spatial-frequency components, then we find 
that Newton's method has a small (poor) radius of convergence, i.e. it may fail to converge 
even when the initial guess is quite close to the true horizon position. In this context there's 
no sharp demarcation between "low" and "high" spatial frequencies, but in practice we use 
the terms to refer to angular Fourier components varying as (say) cos m6 with m ^ 4 and 
m ^ 8 respectively. 

Using a Monte Carlo survey of initial-guess-error Fourier-coefficient space, we find that 
the radius of convergence for Newton's method falls rapidly with increasing spatial frequency, 
approximately as l/m^^^. A simple "line-search" modification of Newton's method roughly 
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doubles the horizon finder's radius of convergence, and shghtly slows the rate of decline with 
spatial frequency. Using a robust nonlinear-algebraic-equations code to solve the discrete 
H(h) = equations would probably give some further improvement, but we doubt that it 
would change the overall trend. 

Using quantitative convergence tests, we demonstrate that the poor high-spatial- 
frequency convergence behavior of Newton's method is not an artifact of insufficient res- 
olution in the finite difference grid. Rather, it appears to be inherent in the (a) strong 
nonlinearity of the continuum H{h) function for high-spatial-frequency components in h. 
We conjecture that H may be similarly nonlinear in its high-spatial-frequency dependence 
on the inverse-metric components. If so, then the presence of high-frequency gravitational 
radiation might well also impair the convergence of Newton's method, and possibly other 
horizon-finding methods as well. Further investigation of this possibility would be very 
interesting. 

Fortunately, if the horizon finder is being used to update a horizon's position at each time 
step of a time evolution, then the previous time step's horizon position probably provides a 
sufficiently good initial guess for Newton's method to converge well. 

Provided it converges, the Newton's-method algorithm for horizon finding is potentially 
very accurate, in practice limited only by the accuracy of the H(h) finite differencing scheme. 
Using 4th order finite differencing, we demonstrate that the error in the numerically com- 
puted horizon position, i.e. the deviation of h from the true continuum horizon position, 
shows the expected 0((A6')^) scaling with grid resolution AO, and is typically ~10~^(10~^) 
for a grid resolution of A9 = ^(fof )• 

Finally, we have argued that considerable further research is needed to develop algorithms 
for finding or recognizing the outermost apparent horizon in a slice. This is an important 
problem for the numerical evolution of multiple-black-hole spacetimes with the black holes 
excluded from the numerical evolution, but so far as we know no reliable algorithms are 
known for it except in spherical symmetry. 
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APPENDIX A: SYMBOLIC-DIFFERENTIATION JACOBIAN COEFFICIENTS 

In this appendix we tabulate all the nonzero symbolic-differentiation Jacobian coefficients 
for ^^^H{h) and its subfunctions. These are used in the 3d.sd.ls and 3d.sd2.2s Jacobian- 
computation methods. All the coefficients are obtained by straightforward, if somewhat 
tedious, linearizations in the manner of Eq. (|2^), starting from the defining equations noted. 

For ^^^H{h), starting from Eqs. and (16), the coefficients are 
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dB 
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dB 



where 
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dA 



d{d,h) 



dB 
d{d,h) 
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d{dji) 

dD 
d{d,h) 

dA 
d{d^yh) 

dB 



+ - m-^a^duh + {d,gn{duh){d,h) 
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- d,g'^ - {d. In ^)g''' 

- 2{K"' - K^^'duh) 
-2{g''^ -g^^duh) 



d{d^yh) 

For Si{h), starting from Eq. (0), the coefficients are 



1 if u = X 
otherwise 



d{d^h) 

For n^{h), starting from Eq. (P), the coefficients are 

dn^ g^^ — g^'^duh)(g^^ — g^^dyh) 

d{d^h) ^ ~D^^ D^^ ■ 

For n^{su), starting from Eq. (^), the coefficients are 

dn' {g''^Sk){g'^^si) 
dsu {g^^SkSiy/"^ {g^^SkSiY/"^ 

For ^^^H{si), starting from Eqs. ([T^) and (|T5l), the coefficients are 
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For ^^^H{n^), starting from Eq. (|^), the coefficients are 
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drv^ 

am 

d{d,ny) 



In + 2K^in' 

J 1 if X = y 
1 otherwise 



(A5c) 
(A5d) 

(A5e) 

(A5f) 

(A5g) 

(A5h) 



(A6a) 
(A6b) 



APPENDIX B: DETAILS OF OUR HORIZON-FINDING CODE 

In this appendix we outhne those details of our horizon-finding code relevant to the 
remainder of this paper. 

Our horizon finder implements all the horizon-function and Jacobian-computation meth- 
ods discussed in this paper, as summarized in table |. It's part of a larger 3 + 1 code under 
development, designed to time-evolve an asymptotically fiat axisymmetric vacuum space- 
time containing a single black hole present in the initial data. The black hole is excluded 
from the numerical grid in the manner described by Refs. [0,^,^,^ . The code uses 4th or- 
der centered finite differencing (5-point molecules) for finite differencing, on a 2-dimensional 
polar-spherical-coordinate grid. (The code also assumes equatorial reflection symmetry, 
but this is merely for convenience and could easily be changed.) The code uses a "PDE 
compiler" to automatically generate all the finite differencing and other grid-computation 
code, including that for the horizon function and Jacobian computations, from a high-level 
tensor-differential-operator specification of the 3 + 1 equations. 

The entire code is freely available on request from the author, and may be modified 
and/or redistributed under the terms of the GNU Public License. The code should be 
easily portable to any modern computing platform. It's mainly written in ANSI C (about 
30K lines) and the Maple symbolic-computation language (about 9K lines for the PDE 
compiler itself, and about 6K hues for the 3 + 1 equations), together with about IK lines 
of awk. The code for the horizon finder itself is about 6K lines of C and 2K lines of 
Maple, but a large part of this is due to its supporting many different combinations of 
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finite differencing schemes and horizon-function and Jacobian computation methods. We 
estimate that an implementation supporting only a single differencing scheme and horizon- 
function and Jacobian computation method, supplemented by a not-optimized-for-efficiency 
independent Jacobian computation for debugging purposes (cf. section [VI D|) , would be a 
factor of ~4 smaller. 

The code takes the metric, extrinsic curvature, and other 3 + 1 field tensors to be al- 
gebraically fully general, i.e. it permits all their coordinate components to be nonzero. To 
avoid z axis coordinate singularities, the code uses a hybrid of polar spherical and Cartesian 
coordinates as a tensor basis. As discussed in detail by Ref. [Q, for the subset of the slice 
containing the code's (2- dimensional) grid, this hybrid coordinate system combines the con- 
venient topology of polar spherical coordinates with the singularity-free nature of Cartesian 
coordinates. 

For present purposes, the key consequence of this 2;-axis-handling method is that in this 
work we've made no effort to avoid expressions which would be singular on the z axis if 
polar spherical coordinates were used as a tensor basis. We haven't investigated this case in 
detail, but we suspect such singularities would be widespread. 



APPENDIX C: OUR SPARSE-JACOBIAN STORAGE SCHEME 



As mentioned in section |VI D| , all the Jacobian matrices involved in horizon finding are 



highly sparse, and for reasonable efficiency this sparsity must be exploited in storing and 
computing the Jacobians. In this appendix we briefly describe our sparse- Jacobian storage 
scheme. This scheme stores the Jacobian by rows, and is applicable to all of the Jacobian 
matrices which arise in our horizon-finding algorithm. 

We consider first the storage of J *^^^H(h) . Which elements in a specified row i of 
this Jacobian are nonzero? From the basic definition Eq. (^), we see that the nonzero 
elements j are precisely those where ^^^Hj depends on hj, i.e. those for which hj enters into the 
computation of *^^)Hi. That is, for a l(2)-stage *^^)H(h) computation, the nonzero- Jacobian j 
values are precisely those within 1(2) molecule radia of i. This makes it easy to store the 
Jacobian: for each grid point i, we simply store a molecule-sized (twice-molecule-sized) array 
of Jacobian elements. 

In practice, for an axisymmetric-spacetime (2-dimensional) code, where i and j are both 
1-dimensional {&) grid indices and the Jacobian is a band matrix, we would store the Jacobian 
as a 2-dimensional array with indices i and j— i. For a fuUy-general-spacetime (3- dimensional) 
code, where i and j are both 2-dimensional [Q and 0) grid indices, we would store the Jacobian 
as a 4-dimensional array with indices i^, i^, — ig, and j,/, — i^, where we temporarily use 
subscripts for coordinate components, and where for pedagogical simplicity we ignore the 
artificial grid boundaries at = {0, vr} and = {0, 27r}. 

A similar storage scheme may be used for more complicated Jacobians. For example, 
consider the storage of J ^^)H(h) . Here i is a 3-dimensional grid point index for '^^^H, while 

J is a 2-dimensional grid point index for h. For a l(2)-stage *^^)H(h) computation, the nonzero 
Jacobian elements in a specified Jacobian row i are now precisely those j within 1(2) angular 
molecule radia of the angular components of i. Thus for an axisymmetric-spacetime (2- 
dimensional) code we would store this Jacobian as a 3-dimensional array with indices i^. 
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i0, and Je — 10, while for a fully-general-spacetime (3-dimensional) code we would store the 
Jacobian as a 5-dimensional array with indices i^, ig, i<^, Je — le, and — i,^. 

Notice that with this storage scheme the Jacobian's structure, i.e. the location of its 
nonzero elements, is stored implicitly. This makes this scheme considerably more efficient 
in both space and time than generic "sparse matrix" storage schemes (for example those 



of Refs. ||3T1 , |39|1 ) , which invariably require the storage of large integer or pointer arrays to 



record a sparse matrix's structure. 



APPENDIX D: DETAILS OF OUR HORIZON-PERTURBATION SURVEY 

In this appendix we describe our Monte Carlo horizon-perturbation survey (cf. sec- 
tion [IX B| ) in more detail. Given the maximum initial-guess-error spatial frequency M, 
the goal of the survey procedure is to estimate Vm, the (hyper) volume in {cm}-space of the 
horizon finder's convergence region. 

To do this, we first start from the origin in {cm}-space, and search outwards along each 
Cm axis until we find coefficients for which the horizon finder fails to converge. This gives 
the intersection of the Cm coordinate axes with the boundary of the convergence region. 

We then construct a sequence of nested hypercubes (strictly speaking, hyper-parallele- 
pipeds) Ci, C2, C3, ... in {cm}-space, starting with Ci just containing the Cm-coordinate- 
axis boundaries of the convergence region, and expanding outwards. We use the obvious 
Monte Carlo sampling algorithm to estimate the volume of the convergence region contained 
within the first hypercube Ci, and then within the differences Ck+i — Ck of the succeeding 
hypercubes. We continue this process until one of the differences contains no convergence- 
region volume. We include one additional hypercube in the sequence after this, typically 
25-50% larger than the previous one in each dimension, to provide a safety margin against 
missing disconnected "islands" or fractal zones near the boundary of the convergence region. 
(These are quite plausible; recall that the (fractal) Julia set is just the convergence region 
of a simple Newton's- method iteration.) Finally, we compute an estimate for Vm by simply 
adding the convergence-region- volume estimates for Ci and each C^+i — Ck- 

Unfortunately, as M and hence the dimensionality of {cm}-space increases, we find that 
the fraction of the hypercubes and hypercube differences occupied by the convergence re- 
gion decreases rapidly, so a very large number of horizon-finding trials is needed to obtain 
a reasonable statistical accuracy for Vm- (For example, the M = 12 points in figure |^ re- 
quired 15,000 trials each.) It's this effect which ultimately limits the maximum value of M 
attainable in practice by a survey of this type. 
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FIGURES 



FIG. 1. This figure illustrates the various 2-stage and 1-stage computation methods for the 
horizon function H(h). The solid arrows denote finite differencing operations, the dotted arrow 
denotes an algebraic computation, and the dashed arrow denotes a radial interpolation to the 
horizon position r = h[6,<j)). Each path from h to H represents a separate computation method. 
Notice that there are three distinct 2-stage methods (using the upper arrows from (2)h to (3)h m 
the figure) and one 1-stage method (using the lower arrow from (2)h to (3)H). 



FIG. 2. This commutative diagram illustrates the two different ways a Jacobian matrix can 
be computed. Given a nonlinear continuum function P{Q), the Jacobian matrix J P(Q) is logi- 
cally defined in terms of the lower-left path in the diagram, i.e. it's defined as the Jacobian of a 
nonlinear discrete (finite difference) approximation P(Q) to P{Q). However, if the operations of 
discretization (finite differencing) and linearization commute, we can instead compute the Jacobian 
by the upper-right path in the diagram, i.e. by first linearizing the continuum P{Q) function, then 
discretizing (finite differencing) this linearization 5P{SQ). 



FIG. 3. This figure shows H{r) for spherical trial horizon surfaces with coordinate radius 
r in an Eddington-Finkelstein slice of a unit-mass Schwarzschild spacetime. Notice that for 
r > r™''^ ~ 4.372, H > and dH/dr < 0, so Newton's method diverges in this region. 



FIG. 4. This figure illustrates how the convergence behavior of the basic and modified Newton 
iterations depends on the spatial-frequency spectrum of the initial guess's error h(o) - h*. In each 
part of the figure, the true continuum horizon h* is plotted as a solid line, while the horizon finder's 
first few iterates (trial horizon surfaces) h^*^^ are plotted with dots at the grid points. Part (a) of 
the figure shows the behavior of Newton's method for an initial-guess-error containing only low 
spatial frequencies, part (b) shows the behavior of Newton's method for an initial-guess-error 
containing significant high spatial frequencies, and part (c) shows the behavior of the modified 
Newton iteration for the same initial guess as part (b). In parts (a) and (c), where the iteration is 
converging, the final iterates shown are indistinguishable from the true continuum horizon at the 
scale of the figure. In part (b), where the iteration is diverging, the computed values for the next 
iterate h^^) (not shown) are almost all far outside the scale of the figure; many of them are in fact 
negative! 



FIG. 5. This figure shows the results of a 3-grid covergence test for the 2nd- iteration Newton 
iterate (trial horizon surface) h^^^ plotted in figure §(b). The line has slope j^, appropriate for 
4th order convergence. (Recall that this line isn't fitted to the data, but is rather an a priori 
prediction with no adjustable parameters.) (The absolute magnitude of the errors shown here is 
much larger than is typical for our horizon finder, due to a combination of the compounding of 
smaller errors in the earlier Newton iterate h'^^^ and the very strong angular variation in both 
iterates h(i) and h(2).) 
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FIG. 6. This figure shows the volume ratios Rm for the horizon-perturbation survey. These 
measure the average radius of convergence of the horizon finder as a function of the ini- 
tial-guess-error's maximum spatial frequency M. The points and solid lines show the results for 
the modified-Newton (upper) and Newton (lower) algorithms, with it la statistical error bars from 
the Monte Carlo estimation procedure. The dashed line shows an Rm ~ falloff. 



FIG. 7. This figure illustrates the accuracy of our horizon finder for a test case where the 
horizon's coordinate shape is strongly non-spherical. The figure is plotted using the transformed 
radial coordinate defined by Eq. (46). Part (a) of the figure shows the "peanut-shaped" true 
continuum horizon position h*, plotted as a solid line, and the initial guess h^'^) and the final 
numerically computed horizon position) h^^'^ plotted with dots at the grid points. At this scale the 
numerically computed horizon position h^^^ is indistinguishable from the true continuum position 
h* . Part (b) of the figure shows the results of a 2-grid convergence test for the numerically computed 
horizon position h^P\ The line has slope jq, appropriate for 4th order convergence. (Recall again 
that this line isn't fitted to the data, but is rather an a priori prediction with no adjustable 
parameters.) 
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TABLES 



TABLE L This table summarizes the various methods for computing the horizon function 
(2)H(h) and its Jacobian J (2)H(h) . The "codes" are shorthand labels for referring to the var- 
ious methods. The relative CPU times are as measured for our implementation (described in 
appendix and are per angular grid point, normalized relative to the 1-stage (^)H(h) computa- 
tion. The notation "Si|n*" means whichever of Sj and/or n* is appropriate, depending on the precise 
2-stage method used to compute the horizon function. 



Code 



H.ls 
H.2s 



2d.np.ls 
2d.np.2s 
3d.np.ls 
3d.sd.ls 
3d.np.2s 
3d.np2.2s 

3d.sd2.2s 



Jacobian 

Computation 

Dimensions 



Jacobian Computation Method 



Horizon 

Function 

Method 



1- stage 

2- stage 



2-dimensional 

2- dimensional 

3- dimensional 
3-dimensional 
3-dimensional 
3-dimensional 

3-dimensional 



Numerical perturbation of (^)H(h) 1-stage 

Numerical perturbation of (^)H(h) 2-stage 

Numerical perturbation of (^)H(h) 1-stage 

Symbolic differentiation of *^^)H(h) 1-stage 

Numerical perturbation of (^)H(h) 2-stage 

Numerical perturbation 2-stage 
of Si|n^(h) and (3)H(si|n') 

Symbolic differentiation 2-stage 
of Si|n*(h) and (3)H(si|nO 



Jacobian Matrices Used 



(2)H(h 
(2)H(h 
(2)H(h 
(2)H(h 
(2)H(h 
(2)H(h 
Si|n*(h 
(2)H(h 



J Si|n*(h 



, j[(3)H(h) 
, j[(3)H(h) 
, j[(3)H(h) 
, j[(3)H(h) 
, j[(3)H(s,|n' ) 
, j[(3)H(h)], 

, jR3)H(s,|nO 
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TABLE I. This table summarizes the various methods for computing the horizon function 
(2)H(h) and its Jacobian j[(2)H(h)]. The "codes" are shorthand labels for referring to the var- 
ious methods. The relative CPU times are as measured for our implementation (described in 
appendix B), and are per angular grid point, normalized relative to the 1-stage (^)H(h) computa- 
tion. The notation "sj|n"' means whichever of Sj- and/or n' is appropriate, depending on the precise 
2-stage method used to compute the horizon function. 





Jacobian 




Horizon 










Relative 


Estimated 




Computation 




Function 










CPU 


Implement ation 


Code 


Dimensions 


Jacobian Computation Method 


Method 


Jacobian Mat 


rices Used 


Time 


Effort 


H.ls 






1-stage 










= 1 


Moderate 


H.2s 






2-stage 










0.7 


Low 


2d.np.ls 


2- dimensional 


Numerical perturbation of (^)H(h) 


1-stage 


J 


"(2)H(h) 






6 


Low 


2d.np.2s 


2- dimensional 


Numerical perturbation of (^)H(h) 


2-stage 


J 


"(2)H(h) 






8 


Low 


3d.np.ls 


3- dimensional 


Numerical perturbation of (■^)H(h) 


1-stage 


J 


"(2)H(h) 


], J 


(3)H(h)] 


7 


Low - Moderate 


3d.sd.ls 


3- dimensional 


Symbolic differentiation of (■^)H(h) 


1-stage 


J 


"(2)H(h) 


], J 


(3)H(h)] 


1.5 


Moderate 
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